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ABSTRACT 



Context. Time-dependent, 3D radiation transfer calculations are important for the modeling of a variety of objects, from supernovae 
and novae to simulations of stellar variability and activity. Furthermore, time-dependent calculations can be used to obtain a 3D 
radiative equilibrium model structure via relaxation in time. 

Aims. We extend our 3D radiative transfer framework to include direct time dependence of the radiation field; i.e., the dl/dt terms are 
fully considered in the solution of radiative transfer problems. 

Methods. We build on the framework that we have described in previous papers in this series and develop a subvoxel method for the 
dl/dt terms. 

Results. We test the implementation by comparing the 3D results to our well tested ID time dependent radiative transfer code in 
spherical symmetry. A simple 3D test model is also presented. 

Conclusions. The 3D time dependent radiative transfer method is now included in our 3D RT framework and in PHOENIX/3D. 
Key words, radiative transfer - methods: numerical - stars: supernovae: general 



1. Introduction 

Supernovae of all types undergo a rather rapid evolution after 
their explosion. During the free-expansion phase, observations 
show fast evolving light curves and changing spectral features. 
In type la supernovae, the radioactive decay of 56 Ni heats the 
envelope, causing a peak in the optical light curve about 20 days 
after the explosion. We have already modeled light curves of 
SNe la with our I D spherically symmetric model atmosp here 
code PHOENIX/ ID (iJack et alJ201 li.l2012l:IWang et alJ2012h . To 
compute more accurate model light curves, the hydrodynamical 
evolution during the free expansion phase needs to be calculated 
in 3D, including the time dependence of the radiation field. This 
requires a time-dependent solution of the 3D radiative transfer 
equation, including the velocity field. In addition, even for the 
calculation of static and stationary 3D atmospheres, time relax- 
ation to radiative equilibrium is a possible method for modelling 
3D atmospheres in energy equilibrium. 

In a seri es of papers presen t ing a 3 D radiative transfer 
frame wo rk (lHauschildt & Baronl 20061 Baron & Hauschildt] 
2007 1 ; lHauschildt & Baronl 120081 120091: iBaronetalJ 
2009; Hausc hildt & Baronl 120101; ISeelmann et alj T2010; 



time dependence in the 3D RT framework using a number of 
test calculations that are discussed below. 

In the following section, we describe the method we use to 
solve the time-dependent 3D radiative transfer. The test calcula- 
tions are presented in section 3 to verify that the implementation 
functions correctly. In the tests, we include scattering and inves- 
tigate atmospheres with time-dependent inner boundary condi- 
tions. 



2. Transfer equation 

IChen et ail d2007h present an approach to solve the radiative 
transfer equation in a flat space time and in the comoving frame. 
The radiative transfer e quation written in t erms of an affine pa- 
rameter Eq. (18) from lChen et alJ (120071) . is given by 



dh (dA\dh ( h 5AA\ h 



(1) 



Hauschild t & Baronll201 ll) . several radiative-transfer test prob- 



lems for several scenarios were considered, always under the 
assumption of time independence. In this paper, we extend our 
framework to include direct time dependence in the solution of 
radiative transfer by considering the dl/dt terms in the 3D radia- 
tive transfer e quation. We used our ID time-dependent radiative 
transfer code dJack et ai]l2009l) to verify the implementation of 



The description of the ^//^/{-discr etization for homolo gous ve- 
locity fields has been presented in iBaron et alj d2009l) . We use 
this discretization method and extend it to include the dl/dt time 
dependence in the solution of the radiative transfer. 

The time-dependent 3D radiative transf er equation along a 
characteristic is a modification of Eq. (15) in lBaron et alJ {2009) 
and given by 



dh, d d 
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where the factor a(t) is simply 



1 



a(t) = -. (3) 
c 

The path length along a characteristic is represented by s, the 
intensity 1^, is a function of the wavelength A and time t along 
the characteristic. F or the detailed derivation of f(s) and a(s) see 
IBaron et all d2009l) . 

The fully implicit discretization of equation ( 1 ) in both wave- 
length and time is given by 



ds 



a(s) Al + ^ + 4a(s) + Xi f(s) 
A/ - Ai-i At 



a(s)- — +T] A f(s). 



(4) 



Following Bar on et al . (2009), this leads to a modification of the 
effective optical length^, which is now defined by 

a(s)Ai a(t) 

+ ~a7 



dr = - \XA.f{s) + 4a(s) + 



ds = ~xds. 



(5) 



A/ - 

IBaron et"al] d2009l) chose an unusual sign convention for x, 
which we have restored to its conventional choice. The modi- 
fied source function S a is then defined by 



dr 



Xa 
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(6) 



(7) 



Our restoration of the conventional sig n choice fo r y alte rs the 
sign of Eq.Qas compared to Eq. 21 of Baro n et al.l d2009t) . 

This approach of including time dependence in the 3D RT 
framework is similar t o the first discre tization method for the 
ID case as described in iJack et al.l d2009l) . The discretization of 
the dl/dt term modifies the generalized optical depth and adds 
an additional term to the generalized source function. Please 
note that the underlying frame differs in the ID and 3D cases. 
In the ID case all momentum quantities are measured by a co- 
moving observer, whereas in the 3D case, only the wavelength 
is measured by a comoving observer, the momentum angles 
are measured in th e observer's frame. Additionally as noted in 
IBaron et"al] ((2009), the path length differs in the 3D case and the 
ID case . In the ID case the path length dsM is that defined by 
MihalU (119801) . which is not a tr ue distance, whereas our ds is 
a true distance (Chen et al. 2007). These differences lead to the 
differences in the equations and in the factor a(t). 

3. Test of the implementation 

The new extension for time dependence in the solution of the 
3D radiative transfer needs to be tested and compared to the 
results of our time-dependent, ID spherically symmetric radia- 
tive transfer code. We implemented the extension as described 
above into our MPI parallelized Fortran 95 code described in 
the previous papers of this series. For all our test calculations, 
we use a sphere with a gray continuum opacity parameterized 
by a power law in the continuum optical depth t s ,j. We inter- 
polate the T st d profile from the ID grid onto the 3D grid using 
two-point power-law interpolation. The opacity on the 3D grid 
is then given by x — — At/At. We use a 3D spherical coordi- 
nate system because the results can be directly comp ared to our 
ID sp herically symmetric radiation transport code (Hauschildt 
1 1 992b . The basic model parameters are 




Fig. 1. A sketch of two inner voxels and the characteristics of 
one particular angle, Iq and represent the resulting specific 
intensities of two characteristics that hit the same voxel. 



1. Inner radius of R c = 5 x 10 10 cm and an outer radius of 
R out = 1 x 10 11 cm, 

2. Optical depth in the range of r„„„ = 10" 6 to T mta = 5, 

3. Grey temperature structure with T model = 10 4 K, 

4. Outer boundary condition 17 = and diffusive inner bound- 
ary condition. 

5. We assume an atmosphere in LTE, 

6. The atmosphere is static. 

For our test calculations, we use a moderately sized 3D grid 
with n r = 2 * 128 + 1, = 2*8 + 1 and ng = 2 * 4 + 1 points 
along each axis, for a total of 257 * 9 * 17 a 2 x 10 4 voxels. 
See section l3~2l for a detailed explanation why this voxel setup 
is used. In all tests, we use the ful l characteristic method fo r the 
3D time-dependent RT solution (Hauschildt & Baronl2006l) . For 
the solid angle sampling (0, (/>), we chose a grid of 64 2 points, 
whic h is a reasonable resolution of the spherical coordinate sys- 
tem (iHauschil dt & Baronll2009l) . The time-dependent test calcu- 
lations are performed with a time step of 10 _2 s, unless stated 
otherwise. 



3. 1 . Constant atmosphere 

First, we verify the results for the case of a model atmosphere 
completely constant in time. The results of the time-dependent 
RT solution should relax in time to be equal to the time- 
independent solution. 

To obtain a numerically accurate solution, it is necessary to 
follow the time derivative of the intensities for each character- 
istic individually (hereafter: subvoxel method). Previously, the 
mean intensity of each voxel is a voxel average of the intensi- 
ties of all characteristics that go through the particular voxel. 
The characteristics are started either from the innermost or out- 
ermost voxels. One specific inner voxel can be hit by many char- 
acteristics of the same angle. Figure Q] shows an example of two 
innermost voxels with the characteristics of one particular angle. 
Clearly, the voxels are hit by the characteristic that started at the 
neighboring voxel. This second hit gives an additional and dif- 
ferent result for specific intensity for the particular voxel. The 
intensities are taken at a point of the characteristic that is closest 
to the center of the voxel. Previously, the mean of all of these 
intensities has been used to compute the intensity of this partic- 
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Fig. 2. Demonstration of numerical errors. The solid line shows 
the mean intensity of one outermost voxel for a method that 
tracks the intensity time derivatives by averaging over a voxel. 
The mean intensity is inaccurate and then relaxes to the wrong 
solution. The dot-dashed line shows the results of the subvoxel 
method, which is numerically much more accurate. 



ular voxel, / mea n = 2 hln. Using the average of the intensities 
passing through a given voxel to model the time-dependent in- 
tensities results in significant numerical errors, and must there- 
fore be avoided. The explanation is that the difference between 
the averaged intensity and the individual intensity of a charac- 
teristic adds an additional term to the source function. This is 
completely analogous to the correct treatment of the wavelength 
derivative in Lagrange frame radiative transfer calculations. In 
the subvoxel method, all the intensities are saved and used for 
the individual characteristics. The subvoxel method increases the 
storage requirements, but this is addressed with a solid-angle do- 
main decomposition method so that the storage requirements per 
process remain limited. 

Figure [2] shows the mean intensity J of one outermost voxel 
as function of time. Although the atmosphere structure is time 
independent, the mean intensity changes in time and relaxes to a 
wrong solution if a numerical method with voxel-averaged time 
derivatives is used. The errors are on the order of 15% compared 
to the time independent solution. The mean intensity obtained 
with the subvoxel method is constant in time and identical to the 
3D time-independent solution. Thus, we only use the subvoxel 
method. 

3.2. Perturbations 

To see direct effects of the time dependence in the solution of 3D 
RT problems, we use a setup with a time variable inner bound- 
ary condition of our test model atmosphere. These perturbations 
will then propagate through the model atmosphere by radiation 
transport in time and finally emerge at the surface of the sphere. 
We perform the same calculations with our ID spherically sym- 
metric code and compare the results to the time-dependent 3D 
spherically symmetric RT. 

For all our perturbation tests, we place a "light bulb" at the 
center of the model atmosphere. This light bulb is just a central 
source of light inside of the atmosphere. For the first test, the 
light bulb is switched on and set to (arbitrarily) radiate at 10 5 
times the initial inner boundary intensity. The calculations then 
track the changes in the 3D (or ID) radiation field in time. 
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Fig. 3. The error of the mean intensity versus t for different num- 
bers of radial voxels, n r , of a model atmosphere with a light bulb 
at the inner boundary that has been switched on at the start of 
the calculations. 



First, we compare the resulting mean intensity J of the time- 
independent solution to the relaxed time-dependent solution, 
which should be identical. In Fig. [3] the error in the resulting 
mean intensity, (Jd ep - Jtndep) / Jindep, is shown after the atmo- 
sphere has relaxed to the new inner boundary condition. Clearly 
the error is lower when a higher resolution in the radial voxels 
n r is used. This is equivalent to a higher resolution in the optical 
depthr. Withn,. = 2*512 + 1, the erroris reduced to below0.1%. 
Therefore, we use a large number of voxels for the radial coor- 
dinate. Changing the number of voxels for the other coordinates 
does not affect the resulting mean intensity in this spherically 
symmetric test. 

The test model run with n r = 2*512+1 has been computed on 
a supercomputer, where we used 2,048 CPUs. The computation 
time for this calculation with 500 time steps is about two hours. 
The storage requirement is about 1 1 GB and is mainly required 
for saving all the intensities at a time step. By using a solid- 
angle domain decomposition, the memory per process required 
is kept small. For true 3D models the resolution in both the other 
coordinates ng and likely also needs to be significantly higher. 
These calculations are beyond the scope of this work. In future 
work, we will test our code with a higher resolution in ng and n$ 
and apply realistic 3D models to further test our time-dependent 
extension. 

In Fig. |U the mean intensity at different layers (radii in the 
3D atmosphere) is shown versus time. The additional radiation 
of the inner light bulb needs time to propagate through the at- 
mosphere. One simple test is to calculate the time the radiation 
would need at the speed of light to move from the inner to the 
outermost layer. For the distance of R ou , - R c as 5 • 10 10 cm, the 
radiation needs about « 1.7 s, which is consistent with the result 
of the 3D radiative transfer as seen in Fig. |4] Figure [5] illustrates 
the propagation of the radiation "wave" throughout the atmo- 
sphere. Here, the mean intensity is plotted against radius for a 
few snapshots at different times. 

Another test is to check if the time relaxed result of the time- 
dependent 3D RT depends on the size of the time steps. The 
results of this test are shown in Fig. [6] the resulting mean inten- 
sities of the outermost layer are plotted for the test case with the 
light bulb inside computed with three different time step sizes. 
The change in the mean intensity emerges at the same point in 
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Fig. 4. The mean intensity versus time for voxels at different op- 
tical depth, t, of a model atmosphere with a light bulb at the 
inner boundary that has been switched on at the start of the cal- 
culations. 
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Fig. 5. Snapshots of the mean intensity of the atmosphere at dif- 
ferent points instants in time. 
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Fig. 6. Mean intensity of the outermost radius versus time com- 
puted with three different time step choices. 



Fig. 7. The mean intensity of the outermost radius versus time 
computed with 3D and ID time-dependent radiative transfer. 
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Fig. 8. The mean intensity at different optical depths, r, with a 
sinusoidally varying light bulb at the center of the sphere. 



time for all sizes of the time step. However, the shape of the 
change is different. The explanation is that with a shorter time 
step, the step function that moves through the atmosphere is bet- 
ter resolved in time. 

Another important check is to compare the results of the 
time-dependent 3D radiative transfer to the results of our ID 
spherically symmetric radiative transfer results. The mean in- 
tensity of the outermost layer of the 3D and the ID RT time- 
dependent calculation is shown in Fig. [7] The test case simulates 
a light bulb at the inner boundary that has been switched on. The 
change in the mean intensity emerges at the same point in time, 
and the final mean intensity is the same. 

The next interesting test case is to place a sinusoidally vary- 
ing light bulb at the center of the sphere and to follow the radi- 
ation field in time. This leads to an atmosphere where the mean 
intensity should vary sinusoidally everywhere. For this test we 
used a sine wave with a full period of 4 s. In Fig. [8] the result- 
ing mean intensities for different radii are plotted versus time. It 
takes about 2 s before the perturbation has moved outwards to 
affect the outermost layers. The sinusoidally varying mean in- 
tensity is then observed in every layer. The phase shift between 
the inner and outer radii is approximately n, which is about 2 s 
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Fig. 9. The mean intensity at different optical depths, t, corrected 
for the radial time delay. The mean intensities are shifted arbi- 
trarily to better illustrate the smoothing. 




Fig. 10. Snapshot of the mean intensity at different phases of the 
sine wave. 



in time. In Fig. |9j we corrected the resulting sine waves at each 
optical depth for the radial time delay. We also shifted the mean 
intensity arbitrarily to directly overlay the sine waves. With this 
correction for the travel time of the light, there is no additional 
phase shift. The smoothing of the sine wave as it moves through 
the atmosphere is also clearly illustrated. A snapshot of the mean 
intensity at different moments in time is shown in Fig. [10] The 
phase difference between the two snapshots is n. 

3.3. Continuum scattering 

In this section, we investigate the effects of continuum scattering 
on the solution of the 3D time dependent radiative transfer. For 
that, we use a model atmosphere with a low optical depth, so that 
we can more easily "see" the light bulb at the center. Therefore, 
we chose for this test: t„„„ = 1CT 2 and r mflA = 5. 

We solved the radiative transfer problem for this model at- 
mosphere both without scattering and for a scattering dominated 
atmosphere with e = 1CT 4 . A visualization of both spheres is 
shown in Fig. [TT] It shows the image for an observer looking 
down one pole of the coordinate system. The lefthand panel 




Fig. 11. Snapshots of the sphere looking down the pole of the 
coordinate system. In the right panel side no scattering is con- 
sidered, while in the left panel the results with scattering are 
shown. The upper row shows a snapshot at time t — s, and the 
lower row is at the time of minimum of the sine wave. 
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Fig. 12. The mean intensity at one outer voxel as a function of 
time for a sinusoidally varying light bulb at the center of the 
sphere. 



displays the solution for an atmosphere where no scattering is 
considered. The outer voxels are dark and not visible to the ob- 
server, and the sphere shows strong limb darkening. The right- 
hand panel shows the results with scattering included in the so- 
lution of the radiative transfer. The outer voxels of the disk are 
significantly brighter as the radiation from the light bulb is scat- 
tered towards the observer, thus the model showing less limb 
darkening. 

We now let the light bulb inside of the sphere vary sinu- 
soidally. The resulting mean intensity versus time is shown in 
Fig. Q~2] for one outermost voxel. The mean intensity is vary- 
ing sinusoidally, as expected. It takes about 2 s for the radiation 
to travel from the light bulb to the surface. The resulting mean 
intensity for the solution of radiative transfer considering scat- 
tering with € = 10~ 4 is shown in Fig. [13] Again a sinusoidal 
variation in intensity is seen at the surface, but it has a smaller 
amplitude than without scattering. In the presence of scattering, 
it also takes more time for the radiation to move through the at- 
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Fig. 13. Illustration of the time-dependent mean intensity at the 
outer boundary for a model with a sinusoidally varying light bulb 
at the center. In this solution, scattering is considered with e = 
1(T 4 . 

mosphere. A visualization is shown at Fig. QT| The lower row 
shows the apparent disk at a minimum of the sine. For the atmo- 
sphere where scattering is included, the intensity of the voxels at 
the outer disk also varies in time. 

3.4. A 3D test model 

All comparisons with our well-tested ID spherically symmet- 
ric code show that the results of the test models for our im- 
plemented time-dependent 3D radiative transfer framework are 
in good agreement. We now want to verify that our 3D time- 
dependent extension also works for test cases of a fully 3D test 
model atmosphere. Unfortunately, we cannot compare these re- 
sults to our ID code, but we can discuss whether the results are 
reasonable. 

We again use a spherically symmetric setup in the opacity 
and radius for the test model atmosphere. For the optical depth, 
we use a range of r„„„ = 10~ 6 to T max = 5. For the inner radius 
we chose R c = 5 x 10 10 cm and for the outer radius R ou , 
I x 10 11 cm. In this 3D test model, no scattering is considered 
(e=l). 

Since we now want to examine 3D effects, we had to in- 
crease our angular resolution. For the following test calculation, 
we used a setup of «,. = 2*64+1, ng - 2*32+1, and = 2*64+1 
voxels. We now place a Gaussian perturbation of the temperature 
and, therefore, the local source function at an off-center posi- 
tion in the sphere. The center of the Gaussian has a value for 
the source function 50 times higher than the value of the source 
function at the center of the sphere and is located 20 voxels away 
from the center at a radius of 5.69 x 10 10 cm. The source function 
is, therefore, given by 



S = S(r) + 50-S(r)-exp{-[(x- 



-xo) 2 +(y-y ) 2 +(z- 



zo) 2 ]/40},(8) 



with a width for the Gaussian perturbation of a 2 = 40. For the 
first time step, the time-independent model without the perturba- 
tion is solved. After the first time step, the perturbation is intro- 
duced and remains constant for the rest of the calculation. This 
means that the resulting mean intensities J will be constant in 
time after the relaxation process. In Fig. [14] the mean intensity / 
of one ring of the outermost voxels around the sphere is shown 
for different times. The perturbation is located inside the sphere 



Fig. 14. The mean intensity, 7, at a ring of outermost voxels 
(8 = 0) at different times for a test case with an off-center per- 
turbation. 



at <p — 7T and = 0. The mean intensity J at the initial time t = Qs 
is also the result for the non-perturbed model. As expected, the 
voxels closer to the perturbation have a higher mean intensity J. 
The intensity then decreases, the farther away the voxel is from 
the perturbation. On the other side of the sphere, there is no ef- 
fect on the mean intensities of the outer voxels. It is also clear 
that it takes more time for the radiation to reach voxels that are 
farther away from the perturbation. 

The calculation for this 3D test model atmosphere run 
needed about 7,000 CPUh. This shows how demanding a fully 
3D time-dependent atmosphere computation is, especially for 
more realistic models in the future. So far, we have only been 
able to test our time-dependent extension with a gray atmo- 
sphere. For realistic models, we need to solve the non-gray radia- 
tive transfer equation. For a type la supernova model atmosphere 
in LTE, on the order of 10,000 wavelength points are needed in 
the computations. The computation time scales roughly linearly 
with the number of wavelength points. 

4. Conclusion 

We have implemented direct time dependence into our 3D radia- 
tive transfer framework. For best numerical accuracy we used a 
sub voxel method for the discretized time derivative. We have 
also shown that it is important to use high resolution in the ra- 
dial direction. To verify the code, we have calculated the solution 
for a number of simple parameterized ID test cases. The inner 
boundary condition was made time dependent to simulate radia- 
tion waves traveling through an atmosphere. All these perturba- 
tions move through the atmosphere as expected. We also com- 
puted a test case with a sinusoidally varying inner light bulbs. We 
compared the results of the 3D time-dependent radiative transfer 
to the results of our well tested ID spherically symmetric radia- 
tive transfer program and found excellent agreement. In a scat- 
tering dominated atmosphere, it takes more time for the radiation 
to move through an atmosphere. We also calculated a fully 3D 
test case. All tests indicate that the new implementations work 
as intended. 
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